clear all
%% ----- Physical Parameters----------

Ux = 1; % Current Velocity in the X direction. Assumed Constant
Uy=  1;  % Current velocity in the Y direction
h = 50; % height of the bed, assumed constant
g=  9.81; % Gravitational Constant


%% ------------ Simulation Parameters % -----------------
Lx=  100;         % Length in the X direction
Ly = 100;         % Length in the Y direction
Nx = 20;         % number of axial nodes in the X direction
Ny=  20;          % Number of Nodes in the Y direction

steps_x = Nx-1;    % number of axial steps
steps_y = Ny-1;
dx = Lx/steps_x;   % axial step spacing, delta x
dy = Ly/steps_y;
dx2 = 1/(dx*dx);
dy2 = 1/(dy*dy);
dt=  1;
beta=0.5;
%% -------------- Block Matrix Elements------------------------%

[S11,P11] = generateHyperbolic(Nx,Ny,beta*0.5*Ux/dx,beta*0.5*Uy/dy,1.0/dt);
[S22,P22] = generateHyperbolic(Nx,Ny,beta*0.5*Ux/dx,beta*0.5*Uy/dy,1.0/dt);
[S12,P12] = generateElliptic(Nx,Ny,beta*h*dx2,beta*h*dy2, -beta*2*h*(dx2+dy2));
%[S12,P12] = generateElliptic(Nx,Ny,-2/15*h*h*h*dx2,-2/15*h*h*h*dy2,2*(2/15*h*h*h*dx+2/15*h*h*h*dx)+h/3 );


I = eye(Nx*Ny);
Z = zeros(Nx*Ny);
S21 = beta*g*I;

%% -- The Matrix
A = [S11 S12; S21 S22];


%% - Preconditioners
PC(:,:,1) = [P11 Z ; Z  P22];
PC(:,:,2) = [P11 Z; S21 P22];
PC(:,:,3) = [P11 S12 ; Z P22];
PC(:,:,4) = [P11 Z ; S21 P22]*[I inv(P11)*S12; Z I];

%condi(1) = cond(A);
eigsMatrix = eig(A);
%spect(1)= max(abs(eigsMatrix));

plot(eigsMatrix,'*')
pause
%% - Compute the Eigenvalue Spectrum
for i=1:4
    i
    A1(:,:,i) = inv(PC(:,:,i))*A;
    [n,n]= size(A1(:,:,i));
    for ii=1:n
        for jj=1:n
            if A1(ii,jj,i) <=1e-6;
                A1(ii,jj,i)=0;
            end
        end
    end
    %     eigs(A1(:,:,i),1,'lm')./eigs(A1(:,:,i), 1,'sm')
    %     pause

   % condi(i+1,1) = cond(A1(:,:,i));
    eigsPC(:,i) = eig(A1(:,:,i));
  %  spect(i+1,1) =  max(abs(eigsPC(:,i)));
    plot(eigsPC(:,i),'*')
    pause
end

spect


%
% plot(eigsMatrix, '+');
% hold on
% plot(eigsPC, '+');


% A = zeros(Nx*Ny , Nx*Ny);
% for ii=1:Nx
%     for jj= 1:Ny
%         index = (jj-1)*Nx + ii;
%         A(index,index) = cc(ii,jj);
%     end
% end
% for ii=2:Nx-1
%     for jj= 2:Ny-1
%         index = (jj-1)*Nx + ii;
%         A(index,index-1)= ww(ii,jj);
%         A(index,index+1) = ee(ii,jj);
%         A(index,index+Nx) = nn(ii,jj);
%         A(index,index-Nx) = ss(ii,jj);
%     end
% end
%
%
%
% for ii=2:Nx-1
%     jj=1;
%     index = (jj-1)*Nx + ii;
%     A(index,index-1)= ww(ii,jj);
%
%     A(index,index+1) = ee(ii,jj);
%     A(index,index+Nx) = nn(ii,jj);
%
%     jj=Ny;
%     index = (jj-1)*Nx + ii;
%     A(index,index-1)= ww(ii,jj);
%     A(index,index+1) = ee(ii,jj);
%     A(index,index-Nx) = ss(ii,jj);
% end
%
% for jj=2:Ny-1
%     ii=1;
%     index = (jj-1)*Nx + ii;
%     A(index,index+1) = ee(ii,jj);
%     A(index,index+Nx) = nn(ii,jj);
%     A(index,index-Nx) = ss(ii,jj);
%     ii=Nx;
%     index = (jj-1)*Nx + ii;
%     A(index,index-1)= ww(ii,jj);
%     A(index,index+Nx) = nn(ii,jj);
%     A(index,index-Nx) = ss(ii,jj);
%
% end






%index = (j-1)*Nx + i;

%ww = index-1. ee = index+1 .. ss = index-Nx  .. nn = index + Nx